*LIS Code for Regression Loops*
*Code for 8 separate runs*
*Jennifer Hook 31 July 2024*

*JOB 1 Runs Regression and Sample Selection, saving the marginal effect and SE for sex and parenthood among groups*
use $mydata/hook/Sample_Unemp_NILF_Selfemp.dta
replace pi11 = . if Sample_all_emp == 0
by dname, sort: egen ndid = count(pi11) if Sample_all_emp == 1
by dname,: egen i = rank(pi11) if Sample_all_emp == 1
gen pctpil_CY =  i/(ndid + 1) *100 if Sample_all_emp == 1
drop i ndid
replace pctpil_CY = . if Sample_all_emp == 0
replace age = age - 25
sum pi11 pctpil_CY
gen dydx_sex = .
gen dyse_sex = .
gen dydx_mother =.
gen dyse_mother =.
gen dydx_sex_HP = .
gen dyse_sex_HP = .
gen dydx_sex_LP = .
gen dyse_sex_LP = .
gen dydx_sex_HC = .
gen dyse_sex_HC = .
gen dydx_sex_LC = .
gen dyse_sex_LC = .
gen dydx_par_WH = .
gen dyse_par_WH = .
gen dydx_par_WL = .
gen dyse_par_WL = .
gen dydx_par_MH = .
gen dyse_par_MH = .
gen dydx_par_ML = .
gen dyse_par_ML = .
gen dydx_ed_WP = .
gen dyse_ed_WP = .
quietly {
levelsof dname  , local(dnames) 
foreach d in `dnames' {
reg pctpil_CY i.Female#i.tertiary#i.parent i.partnered#i.Female c.age##c.age ///
 if dname == "`d'"
est sto m1
*effect of sex/female *
margins if dname == "`d'" , dydx(Female) post
replace dydx_sex = e(b)[1,2] if dname == "`d'" 
replace dyse_sex = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on women*
margins if dname == "`d'" , dydx(parent) at(Female=(1)) post
replace dydx_mother = e(b)[1,2] if dname == "`d'" 
replace dyse_mother = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (1)) post
replace dydx_sex_HP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (1)) post
replace dydx_sex_LP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (0)) post
replace dydx_sex_HC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (0)) post
replace dydx_sex_LC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (1)) post
replace dydx_par_WH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (0)) post
replace dydx_par_WL = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WL = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (1)) post
replace dydx_par_MH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_MH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (0)) post
replace dydx_par_ML = e(b)[1,2] if dname == "`d'" 
replace dyse_par_ML = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of education on mothers*
margins if dname == "`d'" , dydx(tertiary) at(Female=(1) parent = (1)) post
replace dydx_ed_WP = e(b)[1,2] if dname == "`d'" 
replace dyse_ed_WP = sqrt(e(V)[2,2]) if dname == "`d'" 
est drop m1
}
}
sample 5, count by(dname)
tabstat year Sample_all_heck Sample_fulltime_heck Sample_WH_heck Sample_OCC_heck Sample_Public_heck Sample_EXP_heck, by(dname) stat(mean) nototal format(%9.4f)
tabstat dydx_sex dyse_sex dydx_mother dyse_mother dydx_sex_HP dyse_sex_HP dydx_sex_LP dyse_sex_LP dydx_sex_HC dyse_sex_HC dydx_sex_LC dyse_sex_LC, by(dname) stat(mean) nototal format(%9.4f)
tabstat  dydx_par_WH dyse_par_WH dydx_par_WL dyse_par_WL dydx_par_MH dyse_par_MH dydx_par_ML dyse_par_ML dydx_ed_WP dyse_ed_WP, by(dname) stat(mean) nototal format(%9.4f)

*JOB 2 Runs Heckman, saving the marginal effect and SE for sex and parenthood among groups*
use $mydata/hook/Sample_Unemp_NILF_Selfemp.dta
replace pi11 = . if Sample_all_emp == 0
by dname, sort: egen ndid = count(pi11) if Sample_all_emp == 1
by dname,: egen i = rank(pi11) if Sample_all_emp == 1
gen pctpil_CY =  i/(ndid + 1) *100 if Sample_all_emp == 1
drop i ndid
replace pctpil_CY = . if Sample_all_emp == 0
replace age = age - 25
sum pi11 pctpil_CY
gen dydx_sex = .
gen dyse_sex = .
gen dydx_mother =.
gen dyse_mother =.
gen dydx_sex_HP = .
gen dyse_sex_HP = .
gen dydx_sex_LP = .
gen dyse_sex_LP = .
gen dydx_sex_HC = .
gen dyse_sex_HC = .
gen dydx_sex_LC = .
gen dyse_sex_LC = .
gen dydx_par_WH = .
gen dyse_par_WH = .
gen dydx_par_WL = .
gen dyse_par_WL = .
gen dydx_par_MH = .
gen dyse_par_MH = .
gen dydx_par_ML = .
gen dyse_par_ML = .
gen dydx_ed_WP = .
gen dyse_ed_WP = .
quietly {
levelsof dname  , local(dnames) 
foreach d in `dnames' {
 heckman pctpil_CY i.Female#i.tertiary#i.parent i.partnered#i.Female c.age##c.age, ///
    select(emp = i.Female#i.ykid c.transfer c.otherincome i.educ), if dname == "`d'"
est sto m1
*effect of sex/female *
margins if dname == "`d'" , dydx(Female) post
replace dydx_sex = e(b)[1,2] if dname == "`d'" 
replace dyse_sex = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on women*
margins if dname == "`d'" , dydx(parent) at(Female=(1)) post
replace dydx_mother = e(b)[1,2] if dname == "`d'" 
replace dyse_mother = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (1)) post
replace dydx_sex_HP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (1)) post
replace dydx_sex_LP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (0)) post
replace dydx_sex_HC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (0)) post
replace dydx_sex_LC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (1)) post
replace dydx_par_WH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (0)) post
replace dydx_par_WL = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WL = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (1)) post
replace dydx_par_MH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_MH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (0)) post
replace dydx_par_ML = e(b)[1,2] if dname == "`d'" 
replace dyse_par_ML = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of education on mothers*
margins if dname == "`d'" , dydx(tertiary) at(Female=(1) parent = (1)) post
replace dydx_ed_WP = e(b)[1,2] if dname == "`d'" 
replace dyse_ed_WP = sqrt(e(V)[2,2]) if dname == "`d'" 
est drop m1
}
}
sample 5, count by(dname)
tabstat year Sample_all_heck Sample_fulltime_heck Sample_WH_heck Sample_OCC_heck Sample_Public_heck Sample_EXP_heck, by(dname) stat(mean) nototal format(%9.4f)
tabstat dydx_sex dyse_sex dydx_mother dyse_mother dydx_sex_HP dyse_sex_HP dydx_sex_LP dyse_sex_LP dydx_sex_HC dyse_sex_HC dydx_sex_LC dyse_sex_LC, by(dname) stat(mean) nototal format(%9.4f)
tabstat  dydx_par_WH dyse_par_WH dydx_par_WL dyse_par_WL dydx_par_MH dyse_par_MH dydx_par_ML dyse_par_ML dydx_ed_WP dyse_ed_WP, by(dname) stat(mean) nototal format(%9.4f)

*JOB 3 - Add "full-time" to Heckman as micro-level control*
use $mydata/hook/Sample_Unemp_NILF_Selfemp.dta
keep if Sample_fulltime_heck==1
replace pi11 = . if Sample_all_emp == 0
by dname, sort: egen ndid = count(pi11) if Sample_all_emp == 1
by dname,: egen i = rank(pi11) if Sample_all_emp == 1
gen pctpil_CY =  i/(ndid + 1) *100 if Sample_all_emp == 1
drop i ndid
replace pctpil_CY = . if Sample_all_emp == 0
replace age = age - 25
sum pi11 pctpil_CY
gen dydx_sex = .
gen dyse_sex = .
gen dydx_mother =.
gen dyse_mother =.
gen dydx_sex_HP = .
gen dyse_sex_HP = .
gen dydx_sex_LP = .
gen dyse_sex_LP = .
gen dydx_sex_HC = .
gen dyse_sex_HC = .
gen dydx_sex_LC = .
gen dyse_sex_LC = .
gen dydx_par_WH = .
gen dyse_par_WH = .
gen dydx_par_WL = .
gen dyse_par_WL = .
gen dydx_par_MH = .
gen dyse_par_MH = .
gen dydx_par_ML = .
gen dyse_par_ML = .
gen dydx_ed_WP = .
gen dyse_ed_WP = .
quietly {
levelsof dname  , local(dnames) 
foreach d in `dnames' {
 heckman pctpil_CY i.Female#i.tertiary#i.parent i.partnered#i.Female c.age##c.age i.fulltime, ///
    select(emp = i.Female#i.ykid c.transfer c.otherincome i.educ), if dname == "`d'"
est sto m1
*effect of sex/female *
margins if dname == "`d'" , dydx(Female) post
replace dydx_sex = e(b)[1,2] if dname == "`d'" 
replace dyse_sex = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on women*
margins if dname == "`d'" , dydx(parent) at(Female=(1)) post
replace dydx_mother = e(b)[1,2] if dname == "`d'" 
replace dyse_mother = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (1)) post
replace dydx_sex_HP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (1)) post
replace dydx_sex_LP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (0)) post
replace dydx_sex_HC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (0)) post
replace dydx_sex_LC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (1)) post
replace dydx_par_WH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (0)) post
replace dydx_par_WL = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WL = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (1)) post
replace dydx_par_MH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_MH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (0)) post
replace dydx_par_ML = e(b)[1,2] if dname == "`d'" 
replace dyse_par_ML = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of education on mothers*
margins if dname == "`d'" , dydx(tertiary) at(Female=(1) parent = (1)) post
replace dydx_ed_WP = e(b)[1,2] if dname == "`d'" 
replace dyse_ed_WP = sqrt(e(V)[2,2]) if dname == "`d'" 
est drop m1
}
}
sample 5, count by(dname)
tabstat year dydx_sex dyse_sex dydx_mother dyse_mother dydx_sex_HP dyse_sex_HP dydx_sex_LP dyse_sex_LP dydx_sex_HC dyse_sex_HC dydx_sex_LC dyse_sex_LC, by(dname) stat(mean) nototal format(%9.4f)
tabstat  dydx_par_WH dyse_par_WH dydx_par_WL dyse_par_WL dydx_par_MH dyse_par_MH dydx_par_ML dyse_par_ML dydx_ed_WP dyse_ed_WP, by(dname) stat(mean) nototal format(%9.4f)


*JOB 4 - Add "work hours" to Heckman as micro-level control*
use $mydata/hook/Sample_Unemp_NILF_Selfemp.dta
keep if Sample_WH_heck==1
replace pi11 = . if Sample_all_emp == 0
by dname, sort: egen ndid = count(pi11) if Sample_all_emp == 1
by dname,: egen i = rank(pi11) if Sample_all_emp == 1
gen pctpil_CY =  i/(ndid + 1) *100 if Sample_all_emp == 1
drop i ndid
replace pctpil_CY = . if Sample_all_emp == 0
replace age = age - 25
sum pi11 pctpil_CY
gen dydx_sex = .
gen dyse_sex = .
gen dydx_mother =.
gen dyse_mother =.
gen dydx_sex_HP = .
gen dyse_sex_HP = .
gen dydx_sex_LP = .
gen dyse_sex_LP = .
gen dydx_sex_HC = .
gen dyse_sex_HC = .
gen dydx_sex_LC = .
gen dyse_sex_LC = .
gen dydx_par_WH = .
gen dyse_par_WH = .
gen dydx_par_WL = .
gen dyse_par_WL = .
gen dydx_par_MH = .
gen dyse_par_MH = .
gen dydx_par_ML = .
gen dyse_par_ML = .
gen dydx_ed_WP = .
gen dyse_ed_WP = .
quietly {
levelsof dname  , local(dnames) 
foreach d in `dnames' {
 heckman pctpil_CY i.Female#i.tertiary#i.parent i.partnered#i.Female c.age##c.age c.workhours, ///
    select(emp = i.Female#i.ykid c.transfer c.otherincome i.educ), if dname == "`d'"
est sto m1
*effect of sex/female *
margins if dname == "`d'" , dydx(Female) post
replace dydx_sex = e(b)[1,2] if dname == "`d'" 
replace dyse_sex = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on women*
margins if dname == "`d'" , dydx(parent) at(Female=(1)) post
replace dydx_mother = e(b)[1,2] if dname == "`d'" 
replace dyse_mother = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (1)) post
replace dydx_sex_HP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (1)) post
replace dydx_sex_LP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (0)) post
replace dydx_sex_HC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (0)) post
replace dydx_sex_LC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (1)) post
replace dydx_par_WH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (0)) post
replace dydx_par_WL = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WL = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (1)) post
replace dydx_par_MH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_MH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (0)) post
replace dydx_par_ML = e(b)[1,2] if dname == "`d'" 
replace dyse_par_ML = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of education on mothers*
margins if dname == "`d'" , dydx(tertiary) at(Female=(1) parent = (1)) post
replace dydx_ed_WP = e(b)[1,2] if dname == "`d'" 
replace dyse_ed_WP = sqrt(e(V)[2,2]) if dname == "`d'" 
est drop m1
}
}
sample 5, count by(dname)
tabstat year dydx_sex dyse_sex dydx_mother dyse_mother dydx_sex_HP dyse_sex_HP dydx_sex_LP dyse_sex_LP dydx_sex_HC dyse_sex_HC dydx_sex_LC dyse_sex_LC, by(dname) stat(mean) nototal format(%9.4f)
tabstat  dydx_par_WH dyse_par_WH dydx_par_WL dyse_par_WL dydx_par_MH dyse_par_MH dydx_par_ML dyse_par_ML dydx_ed_WP dyse_ed_WP, by(dname) stat(mean) nototal format(%9.4f)

*JOB 5 - Add "occupation" to Heckman as micro-level control*
use $mydata/hook/Sample_Unemp_NILF_Selfemp.dta
keep if Sample_OCC_heck==1
replace pi11 = . if Sample_all_emp == 0
by dname, sort: egen ndid = count(pi11) if Sample_all_emp == 1
by dname,: egen i = rank(pi11) if Sample_all_emp == 1
gen pctpil_CY =  i/(ndid + 1) *100 if Sample_all_emp == 1
drop i ndid
replace pctpil_CY = . if Sample_all_emp == 0
replace age = age - 25
sum pi11 pctpil_CY
gen dydx_sex = .
gen dyse_sex = .
gen dydx_mother =.
gen dyse_mother =.
gen dydx_sex_HP = .
gen dyse_sex_HP = .
gen dydx_sex_LP = .
gen dyse_sex_LP = .
gen dydx_sex_HC = .
gen dyse_sex_HC = .
gen dydx_sex_LC = .
gen dyse_sex_LC = .
gen dydx_par_WH = .
gen dyse_par_WH = .
gen dydx_par_WL = .
gen dyse_par_WL = .
gen dydx_par_MH = .
gen dyse_par_MH = .
gen dydx_par_ML = .
gen dyse_par_ML = .
gen dydx_ed_WP = .
gen dyse_ed_WP = .
quietly {
levelsof dname  , local(dnames) 
foreach d in `dnames' {
 heckman pctpil_CY i.Female#i.tertiary#i.parent i.partnered#i.Female c.age##c.age i.occb1, ///
    select(emp = i.Female#i.ykid c.transfer c.otherincome i.educ), if dname == "`d'"
est sto m1
*effect of sex/female *
margins if dname == "`d'" , dydx(Female) post
replace dydx_sex = e(b)[1,2] if dname == "`d'" 
replace dyse_sex = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on women*
margins if dname == "`d'" , dydx(parent) at(Female=(1)) post
replace dydx_mother = e(b)[1,2] if dname == "`d'" 
replace dyse_mother = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (1)) post
replace dydx_sex_HP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated parents*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (1)) post
replace dydx_sex_LP = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LP = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on highly educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(1) parent = (0)) post
replace dydx_sex_HC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_HC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of sex/female on low educated CHILDLESS*
margins if dname == "`d'" , dydx(Female) at(tertiary=(0) parent = (0)) post
replace dydx_sex_LC = e(b)[1,2] if dname == "`d'" 
replace dyse_sex_LC = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (1)) post
replace dydx_par_WH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated women*
margins if dname == "`d'" , dydx(parent) at(Female=(1) tertiary = (0)) post
replace dydx_par_WL = e(b)[1,2] if dname == "`d'" 
replace dyse_par_WL = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (1)) post
replace dydx_par_MH = e(b)[1,2] if dname == "`d'" 
replace dyse_par_MH = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated MEN*
margins if dname == "`d'" , dydx(parent) at(Female=(0) tertiary = (0)) post
replace dydx_par_ML = e(b)[1,2] if dname == "`d'" 
replace dyse_par_ML = sqrt(e(V)[2,2]) if dname == "`d'" 
est restore m1
*effect of education on mothers*
margins if dname == "`d'" , dydx(tertiary) at(Female=(1) parent = (1)) post
replace dydx_ed_WP = e(b)[1,2] if dname == "`d'" 
replace dyse_ed_WP = sqrt(e(V)[2,2]) if dname == "`d'" 
est drop m1
}
}
sample 5, count by(dname)
tabstat year dydx_sex dyse_sex dydx_mother dyse_mother dydx_sex_HP dyse_sex_HP dydx_sex_LP dyse_sex_LP dydx_sex_HC dyse_sex_HC dydx_sex_LC dyse_sex_LC, by(dname) stat(mean) nototal format(%9.4f)
tabstat  dydx_par_WH dyse_par_WH dydx_par_WL dyse_par_WL dydx_par_MH dyse_par_MH dydx_par_ML dyse_par_ML dydx_ed_WP dyse_ed_WP, by(dname) stat(mean) nototal format(%9.4f)

*repeat PUB, EXP, Combo*
Job 6: keep if Sample_Public_heck==1  i.public1
Job 7: keep if Sample_EXP_heck==1 c.wexptl
Job 8: keep if Sample_OCC_heck==1 & Sample_fulltime_heck==1 i.fulltime i.occb1


*Women Only*
*JOB 2 Runs Heckman, saving the marginal effect and SE for sex and parenthood among groups*
use $mydata/hook/Sample_Unemp_NILF_Selfemp.dta
replace pi11 = . if Sample_all_emp == 0
by dname, sort: egen ndid = count(pi11) if Sample_all_emp == 1
by dname,: egen i = rank(pi11) if Sample_all_emp == 1
gen pctpil_CY =  i/(ndid + 1) *100 if Sample_all_emp == 1
drop i ndid
replace pctpil_CY = . if Sample_all_emp == 0
replace age = age - 25
sum pi11 pctpil_CY
keep if Female==1
recode nchildren (4/18=4)
tab nchildren
gen dydx_mother =.
gen dyse_mother =.
gen dydx_par_WH = .
gen dyse_par_WH = .
gen dydx_par_WL = .
gen dyse_par_WL = .
quietly {
levelsof dname  , local(dnames) 
foreach d in `dnames' {
 heckman pctpil_CY i.tertiary#c.nchildren i.partnered c.age##c.age, ///
    select(emp = i.ykid c.transfer c.otherincome i.educ), if dname == "`d'"
est sto m1
*effect of parenthood on women*
margins if dname == "`d'" , dydx(nchildren) post
replace dydx_mother = e(b)[1,1] if dname == "`d'" 
replace dyse_mother = sqrt(e(V)[1,1]) if dname == "`d'" 
est restore m1
*effect of parenthood on highly educated women*
margins if dname == "`d'" , dydx(nchildren) at(tertiary = (1)) post
replace dydx_par_WH = e(b)[1,1] if dname == "`d'" 
replace dyse_par_WH = sqrt(e(V)[1,1]) if dname == "`d'" 
est restore m1
*effect of parenthood on low educated women*
margins if dname == "`d'" , dydx(nchildren) at(tertiary = (0)) post
replace dydx_par_WL = e(b)[1,1] if dname == "`d'" 
replace dyse_par_WL = sqrt(e(V)[1,1]) if dname == "`d'" 
est restore m1
est drop m1
}
}
sample 5, count by(dname)
tabstat year dydx_mother dyse_mother dydx_par_WH dyse_par_WH dydx_par_WL dyse_par_WL, by(dname) stat(mean) nototal format(%9.4f)
